# library
library(ape) # manual: https://cran.r-project.org/web/packages/ape/ape.pdf
library(nlme) # manual: https://cran.r-project.org/web/packages/nlme/nlme.pdf
library(MASS)
library(car)
library(devtools)
library(tidyverse)
library(broom)
library(ggplot2)
library(ggpubr)

Global effects across species

Correlations (trend only)

piS without log transformation

piS with log transformation - Linearity improves

rho without log transformation

rho with log transformation - Linearity improves

Although Ne is a component of rho, there’s no correlation. Kendall is used

Infering significance of log(piS) and log(rho) correlations across species while accounting for phylo

# load data with single value for omegas per species
spp_averageRates_rec <- read.table("./Data/1_Dataset/20spp_Grapes_PopGenStatistics+alpha+rho.txt",sep="|", header=TRUE)
spp_averageRates_rec$piS_log <- log10(spp_averageRates_rec$piS)
spp_averageRates_rec$rho_log <- log10(spp_averageRates_rec$rho)

spp_averageRates_rec$Species <- gsub("O. novo","O. novo-ulmi subsp. novo-ulmi",spp_averageRates_rec$Species)
spp_averageRates_rec$Species <- gsub("O. americana","O. novo-ulmi subsp. americana",spp_averageRates_rec$Species)

# Load tree
spp_tree_raw <- read.tree(file = "./Data/0_Taxonomy/20spp_phyliptree.phy", text = NULL, tree.names = NULL, skip = 0, comment.char = "", keep.multi = FALSE)

spp_tree_thinned <- read.tree(file = "./Data/0_Taxonomy/6spp_phyliptreeNamed.phy", text = NULL, tree.names = NULL, skip = 0, comment.char = "", keep.multi = FALSE)

# match label
rename_tree <- function(tree_object) {
  tree_object$tip.label <- gsub("'","",tree_object$tip.label)
  tree_object$tip.label <- gsub("Parastagonospora nodorum","P. nodorum",tree_object$tip.label)
  tree_object$tip.label <- gsub("Venturia inaequalis","V. inaequalis",tree_object$tip.label)
  tree_object$tip.label <- gsub("Penicillium biforme","P. biforme",tree_object$tip.label)
  tree_object$tip.label <- gsub("Ophiostoma novo-ulmi subsp. novo-ulmi","O. novo-ulmi subsp. novo-ulmi",tree_object$tip.label)
  tree_object$tip.label <- gsub("Sphaerulina musiva","S. musiva",tree_object$tip.label)
  tree_object$tip.label <- gsub("Zymoseptoria tritici","Z. tritici",tree_object$tip.label)
  tree_object$tip.label <- gsub("Ophiostoma novo-ulmi subsp. americana","O. novo-ulmi subsp. americana",tree_object$tip.label)
  tree_object$tip.label <- gsub("Ophiostoma ulmi","O. ulmi",tree_object$tip.label)
  tree_object$tip.label <- gsub("Cercospora beticola","C. beticola",tree_object$tip.label)
  tree_object$tip.label <- gsub("Botrytis cinerea","B. cinera",tree_object$tip.label)
  tree_object$tip.label <- gsub("Sclerotinia sclerotiorum","S. sclerotiorum",tree_object$tip.label)
  tree_object$tip.label <- gsub("Pyricularia oryzae 4091-5-8","M. oryzae Triticum",tree_object$tip.label)
  tree_object$tip.label <- gsub("Neurospora discreta","N. discreta",tree_object$tip.label)
  tree_object$tip.label <- gsub("Aspergillus flavus","A. flavus",tree_object$tip.label)
  tree_object$tip.label <- gsub("Fusarium graminearum","F. graminearum",tree_object$tip.label)
  tree_object$tip.label <- gsub("Cryphonectria parasitica","C. parasitica",tree_object$tip.label)
  tree_object$tip.label <- gsub("Pyricularia oryzae 70-15","M. oryzae rice",tree_object$tip.label)
  tree_object$tip.label <- gsub("Pyrenophora teres f. teres","P. teres f. teres",tree_object$tip.label)
  tree_object$tip.label <- gsub("Rhynchosporium commune","R. commune",tree_object$tip.label)
  tree_object$tip.label <- gsub("Verticillium dahliae","V. dahliae",tree_object$tip.label)
  return(tree_object)
}

spp_tree <- rename_tree(spp_tree_raw)
spp6_tree <- rename_tree(spp_tree_thinned) # used in the modeling for the six selected species


#spp_summ_gene$Species %in% spp_tree$tip.label
#spp_averageRates_rec$Species %in% spp_tree$tip.label # add into the model..
# Function for convenience:
normtest <- function(model) {
  return(shapiro.test(resid(model)))
}
indeptest <- function(model) {
  return(Box.test(resid(model)[order(fitted(model))], type = "Ljung-Box"))
}

Tree for phylogenetic correction

# check tree. This tree is simply taken from NCBI taxonomy
plot(spp_tree)

Variable: omegaA

Model 1: gls using default restricted maximum likelihood

m_omegaA_rec <- gls(omegaA ~ rho_log + piS_log, spp_averageRates_rec, correlation=corGrafen(1, spp_tree, form = ~Species))
hist(resid(m_omegaA_rec))

plot(resid(m_omegaA_rec)~fitted(m_omegaA_rec))

normtest(m_omegaA_rec)
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(model)
## W = 0.92115, p-value = 0.1043
indeptest(m_omegaA_rec)
## 
##  Box-Ljung test
## 
## data:  resid(model)[order(fitted(model))]
## X-squared = 0.0038497, df = 1, p-value = 0.9505
summary(m_omegaA_rec)
## Generalized least squares fit by REML
##   Model: omegaA ~ rho_log + piS_log 
##   Data: spp_averageRates_rec 
##         AIC       BIC   logLik
##   -40.05173 -35.88567 25.02587
## 
## Correlation Structure: corGrafen
##  Formula: ~Species 
##  Parameter estimate(s):
##         rho 
## 7.01202e-10 
## 
## Coefficients:
##                  Value  Std.Error  t-value p-value
## (Intercept) 0.17776268 0.04400075 4.039992  0.0009
## rho_log     0.02443897 0.01254419 1.948230  0.0681
## piS_log     0.02887547 0.01565994 1.843906  0.0827
## 
##  Correlation: 
##         (Intr) rho_lg
## rho_log  0.708       
## piS_log  0.618 -0.071
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -1.4728082 -0.7230900 -0.1659583  0.9016189  1.3924628 
## 
## Residual standard error: 0.04438237 
## Degrees of freedom: 20 total; 17 residual
# calculate the VIF for each predictor variable in the model
vif(m_omegaA_rec)
##  rho_log  piS_log 
## 1.005034 1.005034
# create vector of VIF values
vif_values <- vif(m_omegaA_rec)

# create horizontal bar chart to display each VIF value
barplot(vif_values, main = "GVIF", horiz = TRUE, col = "steelblue", xlim = c(0,6), cex.axis=1, cex.names=0.4) + abline(v = 3, lwd = 3, lty = 2)

## numeric(0)

Variable: omegaNA

Model 1: gls using default restricted maximum likelihood

m_omegaNA_rec <- gls(omegaNA ~ rho_log + piS_log, spp_averageRates_rec, correlation=corGrafen(1, spp_tree, form = ~Species))
hist(resid(m_omegaNA_rec))

plot(resid(m_omegaNA_rec)~fitted(m_omegaNA_rec))

normtest(m_omegaNA_rec)
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(model)
## W = 0.96114, p-value = 0.5668
indeptest(m_omegaNA_rec)
## 
##  Box-Ljung test
## 
## data:  resid(model)[order(fitted(model))]
## X-squared = 0.18366, df = 1, p-value = 0.6682
summary(m_omegaNA_rec)
## Generalized least squares fit by REML
##   Model: omegaNA ~ rho_log + piS_log 
##   Data: spp_averageRates_rec 
##         AIC      BIC   logLik
##   -28.31197 -24.1459 19.15598
## 
## Correlation Structure: corGrafen
##  Formula: ~Species 
##  Parameter estimate(s):
##       rho 
## 0.1438239 
## 
## Coefficients:
##                   Value  Std.Error   t-value p-value
## (Intercept) -0.07605272 0.06662792 -1.141454  0.2695
## rho_log     -0.03612863 0.01788294 -2.020285  0.0594
## piS_log     -0.05376406 0.02256827 -2.382285  0.0292
## 
##  Correlation: 
##         (Intr) rho_lg
## rho_log 0.706        
## piS_log 0.640  0.001 
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -1.52334794 -0.76725710  0.03974815  0.69854392  1.97466635 
## 
## Residual standard error: 0.0663979 
## Degrees of freedom: 20 total; 17 residual
# Multicollinearity in regression analysis occurs when two or more predictor variables are highly correlated to each other, such that they do not provide unique or independent information in the regression model

# A value between 1 and 5 indicates moderate correlation between a given predictor variable and other predictor variables in the model, but this is often not severe enough to require attention. A value greater than 5 indicates potentially severe correlation.

# calculate the VIF for each predictor variable in the model
vif(m_omegaNA_rec)
##  rho_log  piS_log 
## 1.000001 1.000001
# create vector of VIF values
vif_values <- vif(m_omegaNA_rec)

# create horizontal bar chart to display each VIF value
barplot(vif_values, main = "GVIF", horiz = TRUE, col = "steelblue", xlim = c(0,6), cex.axis=1, cex.names=0.4) + abline(v = 3, lwd = 3, lty = 2)

## numeric(0)

Within genome effect of rho

Visual check for omegaA and omegaNA across different rho categories

spp_summ_rec <- read.table("./Data/1_Dataset/20spp_CatRho_rates+piS+SppRho.txt",sep="|", header=TRUE)
#table(spp_summ_rec$Category_rho)
spp_summ_rec$piS_log <- log10(spp_summ_rec$piS)
spp_summ_rec$MeanSpp_rho_log <- log10(spp_summ_rec$MeanSpp_rho)

# center the rho variables
#spp_summ_rec$Category_rho_centered = spp_summ_rec$Category_rho - spp_summ_rec$MeanSpp_rho
# rho plot
#table(spp_summ_rec$Species)

ggplot(subset(spp_summ_rec, Species %in% c("A. flavus", "B. cinera", "C. beticola","C. parasitica","F. graminearum","M. oryzae rice","M. oryzae Triticum","N. discreta","O. novo-ulmi subsp. americana","O. novo-ulmi subsp. novo-ulmi")), aes(Category_rho, omegaA)) +
  geom_point() +
  facet_wrap(~ Species, ncol=2, scales='free')

ggplot(subset(spp_summ_rec, Species %in% c("O. ulmi", "P. biforme", "P. nodorum","P. teres f. teres","R. commune","S. musiva","S. sclerotiorum","V. dahliae","V. inaequalis","Z. tritici")), aes(Category_rho, omegaA)) +
  geom_point() +
  facet_wrap(~ Species, ncol=2, scales='free')  

ggplot(subset(spp_summ_rec, Species %in% c("A. flavus", "B. cinera", "C. beticola","C. parasitica","F. graminearum","M. oryzae rice","M. oryzae Triticum","N. discreta","O. novo-ulmi subsp. americana","O. novo-ulmi subsp. novo-ulmi")), aes(Category_rho, omegaNA)) +
  geom_point() +
  facet_wrap(~ Species, ncol=2, scales='free')

ggplot(subset(spp_summ_rec, Species %in% c("O. ulmi", "P. biforme", "P. nodorum","P. teres f. teres","R. commune","S. musiva","S. sclerotiorum","V. dahliae","V. inaequalis","Z. tritici")), aes(Category_rho, omegaNA)) +
  geom_point() +
  facet_wrap(~ Species, ncol=2, scales='free')

  • For omegaA there are the three patterns:

    • Species with exponential increase in rho VS omegaA (e.g., S. musiva)

    • Species with almost linear increase in rho VS omegaA (e.g., A. flavus)

    • Species with zero omegaA and rho (as omegaA is 0 thus recombination cannot influence it)

  • So a single model cannot fit all these patterns altogether. Thus to continue with a linear modelling, species were removed if:

    • omegaA = 0

    • no enough rho categories (less than 10)

  • For omegaNA remove species with no enough rho categories (less than 10) and C. parasitica and R. commune (these two species have opposite correlation between omegaNA and rho), and remove “P. teres f. teres”, as it was not possible to calculate omegaA here. The final species dataset is equal to the one for omegaA.

After the filtering desscribed above, these are the species dataset used for omegaA and omegaNA:

#table(spp_summ_rec$Species)
remove_spp_omegaA <- c("B. cinera", "C. parasitica", "F. graminearum", "M. oryzae rice", "M. oryzae Triticum", "N. discreta", "O. novo-ulmi subsp. americana", "O. ulmi", "P. biforme", "P. teres f. teres", "R. commune", "S. sclerotiorum", "V. dahliae", "V. inaequalis")

remove_spp_omegaNA <- c("B. cinera", "C. parasitica", "F. graminearum", "M. oryzae rice", "M. oryzae Triticum", "N. discreta", "O. novo-ulmi subsp. americana", "O. ulmi", "P. biforme", "P. teres f. teres", "R. commune", "S. sclerotiorum", "V. dahliae", "V. inaequalis")

spp_summ_rec_omegaA <- spp_summ_rec[!spp_summ_rec$Species %in% remove_spp_omegaA, ]
spp_summ_rec_omegaNA <- spp_summ_rec[!spp_summ_rec$Species %in% remove_spp_omegaNA, ]


ggplot(spp_summ_rec_omegaA, aes(Category_rho, omegaA)) +
  geom_point() +
  facet_wrap(~ Species, ncol=2, scales='free') + stat_cor(method = "kendall") + theme(strip.text = element_text(face = "italic"))

ggplot(spp_summ_rec_omegaNA, aes(Category_rho, omegaNA)) +
  geom_point() +
  facet_wrap(~ Species, ncol=2, scales='free') + stat_cor(method = "kendall") + theme(strip.text = element_text(face = "italic"))

Variable: omegaA - Rho centering method - log

OBS: piS_log is the single value per species

# center the rho variables - method 2
spp_summ_rec_omegaA$Category_rho_log = log10(spp_summ_rec_omegaA$Category_rho)
spp_summ_rec_omegaA$MeanSpp_rho_log = log10(spp_summ_rec_omegaA$MeanSpp_rho)
spp_summ_rec_omegaA$Category_rho_centered = spp_summ_rec_omegaA$Category_rho_log - spp_summ_rec_omegaA$MeanSpp_rho_log

spp_summ_rec_omegaNA$Category_rho_log = log10(spp_summ_rec_omegaNA$Category_rho)
spp_summ_rec_omegaNA$MeanSpp_rho_log = log10(spp_summ_rec_omegaNA$MeanSpp_rho)
spp_summ_rec_omegaNA$Category_rho_centered = spp_summ_rec_omegaNA$Category_rho_log - spp_summ_rec_omegaNA$MeanSpp_rho_log

Model 1: gls using default restricted maximum likelihood

m_omegaA_rec <- gls(omegaA ~ Category_rho_centered + piS_log, spp_summ_rec_omegaA, correlation=corGrafen(1, spp6_tree, form = ~Species))
hist(resid(m_omegaA_rec))

plot(resid(m_omegaA_rec)~fitted(m_omegaA_rec))

normtest(m_omegaA_rec)
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(model)
## W = 0.94368, p-value = 0.001525
indeptest(m_omegaA_rec)
## 
##  Box-Ljung test
## 
## data:  resid(model)[order(fitted(model))]
## X-squared = 4.8557, df = 1, p-value = 0.02756
summary(m_omegaA_rec)
## Generalized least squares fit by REML
##   Model: omegaA ~ Category_rho_centered + piS_log 
##   Data: spp_summ_rec_omegaA 
##         AIC      BIC   logLik
##   -324.1511 -312.432 167.0755
## 
## Correlation Structure: corGrafen
##  Formula: ~Species 
##  Parameter estimate(s):
##          rho 
## 1.238793e-08 
## 
## Coefficients:
##                            Value  Std.Error  t-value p-value
## (Intercept)           0.11713340 0.03431065 3.413908  0.0010
## Category_rho_centered 0.02874998 0.00394184 7.293535  0.0000
## piS_log               0.01191653 0.02027200 0.587832  0.5584
## 
##  Correlation: 
##                       (Intr) Ctgr__
## Category_rho_centered -0.041       
## piS_log                0.949 -0.039
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -3.4701194 -1.7334627 -0.5944074  0.9519679  2.5730211 
## 
## Residual standard error: 0.02655143 
## Degrees of freedom: 80 total; 77 residual
# calculate the VIF for each predictor variable in the model
vif(m_omegaA_rec)
## Category_rho_centered               piS_log 
##              1.001509              1.001509
# create vector of VIF values
vif_values <- vif(m_omegaA_rec)

# create horizontal bar chart to display each VIF value
barplot(vif_values, main = "GVIF", horiz = TRUE, col = "steelblue", xlim = c(0,6), cex.axis=1, cex.names=0.4) + abline(v = 3, lwd = 3, lty = 2)

## numeric(0)

Model 2: gls using default restricted maximum likelihood and bcnPower transformation to improve normality

summary(p1 <- powerTransform(omegaA ~ Category_rho_centered + piS_log, spp_summ_rec_omegaA, family = "bcnPower"))
spp_summ_rec_omegaA$omegaA_trans <- bcnPower(spp_summ_rec_omegaA$omegaA, lambda = p1$lambda, gamma = p1$gamma)
m_omegaA_rec.trans <- gls(omegaA_trans ~ Category_rho_centered + piS_log, spp_summ_rec_omegaA, correlation=corGrafen(1, spp6_tree, form = ~Species))
hist(resid(m_omegaA_rec.trans))

plot(resid(m_omegaA_rec.trans)~fitted(m_omegaA_rec.trans))

normtest(m_omegaA_rec.trans)
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(model)
## W = 0.93155, p-value = 0.0003494
indeptest(m_omegaA_rec.trans)
## 
##  Box-Ljung test
## 
## data:  resid(model)[order(fitted(model))]
## X-squared = 8.0099, df = 1, p-value = 0.004652
summary(m_omegaA_rec.trans)
## Generalized least squares fit by REML
##   Model: omegaA_trans ~ Category_rho_centered + piS_log 
##   Data: spp_summ_rec_omegaA 
##         AIC       BIC   logLik
##   -351.6037 -339.8847 180.8019
## 
## Correlation Structure: corGrafen
##  Formula: ~Species 
##  Parameter estimate(s):
##          rho 
## 1.982586e-08 
## 
## Coefficients:
##                            Value   Std.Error    t-value p-value
## (Intercept)           -0.9040995 0.028708454 -31.492449  0.0000
## Category_rho_centered  0.0236350 0.003298225   7.165963  0.0000
## piS_log                0.0125861 0.016962017   0.742017  0.4603
## 
##  Correlation: 
##                       (Intr) Ctgr__
## Category_rho_centered -0.041       
## piS_log                0.949 -0.039
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -3.5175500 -1.9006669 -0.7021306  1.1196993  2.5516517 
## 
## Residual standard error: 0.02221615 
## Degrees of freedom: 80 total; 77 residual

Variable: omegaNA

Model 1: gls using default restricted maximum likelihood

m_omegaNA_rec <- gls(omegaNA ~ Category_rho_centered + piS_log, spp_summ_rec_omegaNA, correlation=corGrafen(1, spp6_tree, form = ~Species))
hist(resid(m_omegaNA_rec))

plot(resid(m_omegaNA_rec)~fitted(m_omegaNA_rec))

normtest(m_omegaNA_rec)
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(model)
## W = 0.95115, p-value = 0.004012
indeptest(m_omegaNA_rec)
## 
##  Box-Ljung test
## 
## data:  resid(model)[order(fitted(model))]
## X-squared = 0.015695, df = 1, p-value = 0.9003
summary(m_omegaNA_rec)
## Generalized least squares fit by REML
##   Model: omegaNA ~ Category_rho_centered + piS_log 
##   Data: spp_summ_rec_omegaNA 
##         AIC       BIC   logLik
##   -339.9237 -328.2046 174.9618
## 
## Correlation Structure: corGrafen
##  Formula: ~Species 
##  Parameter estimate(s):
##       rho 
## 0.4918982 
## 
## Coefficients:
##                             Value  Std.Error   t-value p-value
## (Intercept)            0.06685267 0.03193614  2.093323  0.0396
## Category_rho_centered -0.02298268 0.00348763 -6.589776  0.0000
## piS_log                0.00783538 0.01714736  0.456944  0.6490
## 
##  Correlation: 
##                       (Intr) Ctgr__
## Category_rho_centered -0.067       
## piS_log                0.925 -0.078
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -2.46673416 -0.70949118 -0.08111046  0.80301592  3.62625371 
## 
## Residual standard error: 0.02400525 
## Degrees of freedom: 80 total; 77 residual
# Multicollinearity in regression analysis occurs when two or more predictor variables are highly correlated to each other, such that they do not provide unique or independent information in the regression model

# A value between 1 and 5 indicates moderate correlation between a given predictor variable and other predictor variables in the model, but this is often not severe enough to require attention. A value greater than 5 indicates potentially severe correlation.

# calculate the VIF for each predictor variable in the model
vif(m_omegaNA_rec)
## Category_rho_centered               piS_log 
##              1.006161              1.006161
# create vector of VIF values
vif_values <- vif(m_omegaNA_rec)

# create horizontal bar chart to display each VIF value
barplot(vif_values, main = "GVIF", horiz = TRUE, col = "steelblue", xlim = c(0,6), cex.axis=1, cex.names=0.4) + abline(v = 3, lwd = 3, lty = 2)

## numeric(0)

Model 2: gls using default restricted maximum likelihood and bcnPower transformation to improve normality

summary(p1 <- powerTransform(omegaNA ~ Category_rho_centered + piS_log, spp_summ_rec_omegaNA, family = "bcnPower"))
spp_summ_rec_omegaNA$omegaNA_trans <- bcnPower(spp_summ_rec_omegaNA$omegaNA, lambda = p1$lambda, gamma = p1$gamma)
m_omegaNA_rec.trans <- gls(omegaNA_trans ~ Category_rho_centered + piS_log, spp_summ_rec_omegaNA, correlation=corGrafen(1, spp6_tree, form = ~Species))
hist(resid(m_omegaNA_rec.trans))

plot(resid(m_omegaNA_rec.trans)~fitted(m_omegaNA_rec.trans))

normtest(m_omegaNA_rec.trans)
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(model)
## W = 0.98345, p-value = 0.3898
indeptest(m_omegaNA_rec.trans)
## 
##  Box-Ljung test
## 
## data:  resid(model)[order(fitted(model))]
## X-squared = 0.0006876, df = 1, p-value = 0.9791
summary(m_omegaNA_rec.trans)
## Generalized least squares fit by REML
##   Model: omegaNA_trans ~ Category_rho_centered + piS_log 
##   Data: spp_summ_rec_omegaNA 
##        AIC      BIC    logLik
##   454.6723 466.3914 -222.3362
## 
## Correlation Structure: corGrafen
##  Formula: ~Species 
##  Parameter estimate(s):
##       rho 
## 0.8026491 
## 
## Coefficients:
##                            Value Std.Error   t-value p-value
## (Intercept)           -18.008024  5.522972 -3.260567  0.0017
## Category_rho_centered  -3.695706  0.599633 -6.163278  0.0000
## piS_log                 1.981752  2.855455  0.694023  0.4898
## 
##  Correlation: 
##                       (Intr) Ctgr__
## Category_rho_centered -0.052       
## piS_log                0.915 -0.077
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.7041156 -0.8992904  0.3256301  1.2707378  3.3866414 
## 
## Residual standard error: 4.179405 
## Degrees of freedom: 80 total; 77 residual
# Multicollinearity in regression analysis occurs when two or more predictor variables are highly correlated to each other, such that they do not provide unique or independent information in the regression model

# A value between 1 and 5 indicates moderate correlation between a given predictor variable and other predictor variables in the model, but this is often not severe enough to require attention. A value greater than 5 indicates potentially severe correlation.

# calculate the VIF for each predictor variable in the model
vif(m_omegaNA_rec.trans)
## Category_rho_centered               piS_log 
##              1.005969              1.005969
# create vector of VIF values
vif_values <- vif(m_omegaNA_rec.trans)

# create horizontal bar chart to display each VIF value
barplot(vif_values, main = "GVIF", horiz = TRUE, col = "steelblue", xlim = c(0,6), cex.axis=1, cex.names=0.4) + abline(v = 3, lwd = 3, lty = 2)

## numeric(0)